/*****************************************************************************/
/*                                                                           */
/*    Domino Operation System Mathematics Module                             */
/*                                                                           */
/*    Copyright (C) 2005 Laszlo Arvai                                        */
/*                                                                           */
/*    ------------------------------------------------------------------     */
/*    mathComplex.h - Complex number calculation routines                    */
/*****************************************************************************/

#include <math.h>
#include "mathComplex.h"

///////////////////////////////////////////////////////////////////////////////
// Types
typedef struct { dosWord low; dosWord high; } THighLow;
typedef union  { dosInt32 lword; THighLow hlword; } TSignedLong;
typedef union  { dosUInt32 lword; THighLow hlword; } TUnsignedLong;

///////////////////////////////////////////////////////////////////////////////
// Local Function prototypes
static dosInt32 IntMult(dosInt32 in_poi1, dosUInt32 in_poi2);

///////////////////////////////////////////////////////////////////////////////
// Complex addition
void mathCplxAdd( dosComplex* in_op1, dosComplex* in_op2, dosComplex* out_result)
{
  out_result->re = in_op1->re + in_op2->re;
  out_result->im = in_op1->im + in_op2->im;
}

///////////////////////////////////////////////////////////////////////////////
// Complex substraction
void mathCplxSub( dosComplex* in_op1, dosComplex* in_op2, dosComplex* out_result)
{
  out_result->re = in_op1->re - in_op2->re;
  out_result->im = in_op1->im - in_op2->im;
}

///////////////////////////////////////////////////////////////////////////////
// Complex multiply for FFT
void mathCplxMultFFT( dosComplex* in_op1, dosComplex* in_op2, dosComplex* out_result)
{
	dosUInt32 cxre,cxim;
  dosInt32 re1,re2,im1,im2;

  cxre = labs(in_op2->re)*2;
  cxim = labs(in_op2->im)*2;

  re1 = IntMult(in_op1->re, cxre );
  re2 = IntMult(in_op1->im, cxim );
  im1 = IntMult(in_op1->re, cxim );
  im2 = IntMult(in_op1->im, cxre );

  if(in_op2->re < 0)
	{ 
		re1 =- re1; 
		im2 =- im2; 
	}

  if(in_op2->im<0)
	{
		re2 =- re2;
		im1 =- im1;
	}

  out_result->re = re1 - re2;
  out_result->im = im1 + im2;
}

///////////////////////////////////////////////////////////////////////////////
// Integer multiply
dosInt32 IntMult(dosInt32 in_poi1, dosUInt32 in_poi2) 
{ 
	dosInt32 l;
  TSignedLong op1,s,t;
	TUnsignedLong poi2;

	poi2.lword = in_poi2;

  op1.lword = labs( in_poi1 );

  s.lword = (dosInt32)op1.hlword.high * poi2.hlword.low;
  t.lword = (dosInt32)op1.hlword.low * poi2.hlword.high;

  l = (dosInt32)op1.hlword.high * poi2.hlword.high +
			(dosInt32)s.hlword.high + (dosInt32)t.hlword.high;

  return( in_poi1 < 0 ? -l : l );
}
